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CO ' Abstract 

<n : 

^3 We study the ground-state properties of dilute Bose gas confined to both isotropic 

^ | and anisotropic traps to assess the accuracy of Gross-Pitaevskii (GP) theory. To 

f— >, ■ go beyond GP approximation we use Huang- Yang theory of interatomic interaction 

^^ , energy for hard-sphere Bose gas and use a variational method to solve the result- 



ed 

S 



ing modified GP equation. We also make an analytic estimate of the corrections 
due to the higher-order terms in the interatomic interaction energy. We find that 
corrections are of the order of 1%. However, there is a qualitative change in the 
density profile due to presence of logarthimic term in the interaction energy for 
large N (number of atoms). 
O 

iri ■ 1 Introduction 

5— i 

The discovery of Bose-Einstein condensation (BEC) in magnetically trapped alkali atoms 

|I], 0, 0] and spin-polarized hydrogen atom || has resulted in large body of theoretical 

and experimental research || on the ground-state as well as excited state properties of 

trapped atomic gases. For the case of dilute atomic gas, which is the case for most of 

the experimental situation, the mean-field theory of GP || has been quite successful 

in explaining both static (ground-state) and dynamic (collective excitations) properties 

of the condensates. The interacting atomic gas is considered to be dilute when the 

average distance between the atoms is much larger than the range of interaction potential 

between the atoms J7| (mathematically this is expressed by the condition n\a\ 3 << 1, 

where n atomic density and a is the s wave scattering length of interatomic potential). 

For this case the interaction between the atoms is well described by two-body collisions 

characterized by a single parameter obtainable from the potential namely, the s wave 

scattering length a. 

In some recent experiments [HI §, [J BEC is achieved with number of atoms N being 

quite high. It is possible to produce condensates of ~ 10 7 87 Rb atoms, ~ 5 x 10 7 Na atoms 

and ~ 10 9 H atoms. The maximum density in the optical traps has been around 10 15 

1 



cm -3 . Thus, with the increasing atomic density in the condensate it becomes necessary 
to test the validity of mean-field GP theory. The main aim of this paper is to investigate 
systematically the accuracy of the GP theory in the high density limit. We go beyond 
mean-field theory and evaluate the corrections using a variational formalism. We also 
derive an analytic estimate for the corrections in the frame work of Thomas-Fermi(TF) 



approximation. For this purpose we use Huang- Yang [JTOJ] theory for the interatomic 
interaction energy e. In this theory, e is expanded in power of gas parameter n\a\ 3 by as- 
suming the atomic gas to be a uniform hard-sphere Bose gas. It is important to note here 
that the interatomic potential in Huang- Yang theory is replaced by a pseudopotential. 
The expression for e is [[LJ], |T^1 

2ixh 2 an 

e[n) = 

m 
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In GP theory only the first term in the above expansion is considered as higher order 
terms are negligible for dilute gas. However with increasing density higher-order terms 
become appreciable. It is, therefore, of interest to study the effects these higher order 
terms on the ground-state properties of BEG 

The effect of higher-order terms in the interatomic potential on the physical properties 
of BEC has recently been reported in the literature [13|, [H]]. Nunes [|13] has used density 
functional theory to derive the modified GP (MGP) equation. By solving the MGP 
equation for a weak isotropic trap in the Thomas-Fermi (TF) approximaion, he finds 
that the maximum change in the ground state density due to the higher order terms in 
Eq.(|l])is only a few percent for the value of N as large as 10 9 . On the other hand, the 
change in the ground state energy per particle has been reported to be sensitive to N - 
for iV < 10 7 both GP and MGP results are very close whereas, they differ by 33% for 
iV = 10 9 . This result is not in conformity with the analytic estimate || of the correction 
introduced in the ground state energy by the higher order terms in Eq. ([!]). According 
to the analytic estimate, the first higher order term in Eq.(|TJ) causes a change of only 
a few percent in the ground state energy even for very large value of N. As discussed 
below, the second higher order term in Eq. ([!]) being logarithmic gives negative correction. 
Therefore, this is expected to partially cancel the correction due to the first higher-order 
term in the interatomic interaction energy. As a result of this, the net correction remain 
close to a few percent when N is varied from 10 7 to 10 9 . Fabrocini and Polls [0 have 



studied the effect of interatomic correlations within the local density approximation and 
the correlated basis function (CBF) theory. They have solved the resulting MGP equation 
and the lowest order correlated Hartree (CH^o) equation for an isotropic trap using the 
steepest descent method. Both of these methods give corrections which are close to those 
predicted by the analytical estimates. However, the CH^o equation gives lower value of 
the ground state energy than the MGP equation. The difference between the results 
obtained by these two methods shows up more at higher values of N . It is not clear if 
this is due to any numerical problem or due to some physical reason. We discuss more 
about these results in section 4. 



In this paper we employ variational approach |I5,|TE, I7|| to solve the MGP equation by 



using the variational ansatz for the density recently proposed in Ref . []T7j . We also perform 
calculations for more realistic case of axially symmetric trap (u x = u y = uj±_ ^ u z ) 
potential. In order to assess the accuracy of our variational results we generalize the 
virial relations || for the interatomic interaction energy given by Eq.([l]). 

Recently, several authors ||, [18], have estimated the corrections introduced by 
the first higher-order term in the Eq.([l]) to assess the accuracy of GP theory. These 
authors obtained analytical expressions for the corrections over mean-field results within 
the Thomas- Fermi (TF) approximation valid in the large N limit. Here we extend these 
results by including second higher-order term of Eq. ([!]). Beside, providing the estimates 
for corrections these results also serve as a check for the accuracy of our numbers obtained 
by the variational calculation. 

It is necessary to emphasize here that our previous and the present study clearly bring 
out the advantages of variational approach. This approach requires numerically less effort 
than the direct numerical integration method and can handle system of arbitrary number 
of particle. At the same time, it provides a very accurate description of the ground state 
properties of the system. 

Rest of the paper is organized as follows: In section 2 we describe the variational 
approach and also derive an expression for the generalized virial relation. Section 3 
contains analytical estimates for corrections in the density profile, the total energy and 
the chemical potential. Results are discussed in section 4 and paper is concluded in 
section 5. 



2 Variational method 

The energy functional for a system of N-bosons, each of mass m, confined in a trap 
potential V t (r) is given by 



E[V] = ( dx 



2m 



2 
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" I W(r)| 2 + V t (vm(v)\ 2 + e(n)|^(r)| 2 



where \l/(r) is the corresponding wave function and gives density via the relation n(r) = 
^(r)! 2 . For e(n) we use the local-density approximation (LDA) expression as given in 
Eq. ([[]). The energy functional associated with GP theory can be obtained from Eq.(|2|) 
by substituting only the first term of Eq.(|l|) for e(n) in Eq.@. In this study we go 
beyond GP theory by including the second and the third term of Eq.(|T]) to study the 
ground-state properties of BEC. 

For our purpose we consider an axially symmetric trap characterized by two angular 
frequencies L/j\ and uo° z and the corresponding potential is given by 

V t (r) = ^f^ + y* + \y), (3) 



where Ao = Wj_/w® is the anisotropy parameter of the trapping potential. 

As in the GP approach, from Eqs.(|]) and (|2|) the MGP equation can be derived by 
making the energy functional E[^f] stationary with respect to the wave function ^(r) 
such that the wave function remains normalized to the total number of particles N given 
by the condition 

"|^(r)| 2 dr = iV. (4) 
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Such a variation leads to: 
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= 2^Vi(ri) (5) 

where, [i\ is the chemical potential in the units of %uf\_ corresponding to the normalization 
constraint given by Eq.(£|). In this equation we have used the length scale rescaled 

variables ri = r/aho,a = a/a ho , Vi = a ho V, E\ = E/hu^_, with a ho = {Ti/muj L )' 1 and 
\I/(ri) is normalized to unity 

In Ref.|n| and ]13[ above equation (Eq.^|) has been solved numerically for isotropic 



trap (Ao = 1)- In this paper we take recourse to variational approach to solve the Eq.(||). 
The main advantage of this method is that with a suitable choice for the variational form 
of the wave function one can get quite accurate results with less computational effort. In 
addition, this approach may provide physical insight which generally get obscured by the 
complicated computational procedures. In this paper we use the variational form given 
by0 




''-('■-) = , ^-T-A-(^) , AW ' T - + X ' T) (6) 



where A, uo± and p are the variational parameters which are obtained by minimizing E\ 

with respect to this parameters. This variational form has been shown to describe quite 

accurately the ground state of the dilute Bose gas confined in a trap for a wide range 

of particle numbers. When number is small it tends correctly towards Gaussian and in 

the opposite limit it resembles with the TF wave function. In addition to this in the 

intermediate region it combines the feature of both in an effective way. This variational 

form is also well suited for negative scattering length and is easily generalized for the 

vortex states. Further, using Eq. @ the physical observables can be expressed analytically 

in terms of these three variational parameters. For example, energy functional E\ can be 

written as: 

E\ 

— = E kin + E ho + E int + E 2 int + E int . (7) 

Here Ekin, Eho denote average kinetic and harmonic trapping potential energies, re- 
spectively. Remaining three terms E} nt , Ef nt and Ef nt represent interatomic interaction 



energies corresponding to the first, second and the third terms in the energy expansion 
inEq.®. 
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For particular value N the parameter uj±, X and p are determined by minimizing the 
energy given above (Eq.(|7|)). 

The virial relation among different energy components provides a way of checking 
the correctness of variational or numerical solutions. For example, within GP theory 
different energy components satisfy the virial relation || 

2E kin — 2E ho + 3E int = (13) 

For our purpose we generalize the above relation (Eq. flllf )) valid for the MGP theory. By 
using variational nature of energy and the scaling transformation ^(r^) — > vAf(r) we 
arrive at following expression for the virial relation corresponding to the MGP energy 
functional 

+ 3EL + 1&L. + 2Ef nt + J |^i( ri )| 6 rf ri = 0. (14) 

For each term of this expansion we have checked the satisfaction of virial relation corre- 
sponding to our variational solutions. The virial relation is satisfied up to 5-th decimal 
place in our calculation. 

Before presenting the results, we make an estimate of the corrections in the density, 
total energy and the chemical potential as a result of inclusion of two new terms in the 
interatomic energy. For this we follow the approach of Dalfovo et al.||, which is valid 
for large N. 



•-kin ~ ±n ho -r oEj int + -E? nt + 2E; 
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3 Analytic Estimate 

We use local density approximation for the chemical potential which is justified in the 
thermodynamic limit N — > oo where the profile of the density distribution is very smooth. 
In this approximation the chemical potential is given by: 



/i = }J>iocai[n{r)] +V t (r) 



(15) 



Given the expression for fii oca i( n ) for the homogeneous Bose-gas, the density n(r) can 
be calculated by solving Eg. fllS]) , along with the normalization condition which fixes the 
parameter \x on the left hand side of the equation. From Eq. ([!]), it follows that 



Hioccd(n) = gn 
where, 
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Considering only the first term in the expansion for fii oca i we recover the mean- field 
Thomas- Fermi result 

n ( r ) = it ~ W r )) 
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n TF {r). 



Normalization condition implies that 



/' 



(15 NX a 



Htf- 



(19) 



Using the relation \x = dE/dN we get the corresponding expression for the ground sate 
energy as 



E = -N/i = E TF . 



(20) 



When the higher order terms in the expansion for \ii oca i are considered the standard 
procedure to solve Eq . fllBD is by iteration ||. On following such procedure we get 

(ji-Vtir)) 4m 3 / 2 
n(r) = : — (fi- V t {r)) ' 



g 3n 2 h 3 



ma 



(fi-V t (v)Y\n l—pfr-Vtiv)) \+0(na d ). (21) 



Now imposing the normalization condition and neglecting the terms of order (9(na 3 )and 
higher, we obtain 



H = Htf 



1 + yW^O) + |j (y ~ Vs\ (n(0)a 3 ) In (n(0)a 3 ) 
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and consequently 

1 + -yVaMO) + ^ ( y " v^J (n(0)a 3 ) In (n(0)a 3 ) 
The parameter a 3 n(0) appearing in the above equations can be written as 



5 
E = -N^tf 



(23) 



1 K~ \ — 

« 8 »(°) = yt r A »«) 5 ( 24 ) 

Thus by using Eqs.(|n]),(JZ2"D and Q2lf ) we estimate the correction over the corresponding 
TF results. For example a BEC with N = 10 7 particles and a = 4.33 x 1CT 3 , the second 
terms in Eq.(|2"2"D and (^) result in 2.2% and 1.9% corrections in the chemical poten- 
tial and the total energy, respectively, over the corresponding TF values. Third terms, 
however, lead to negative corrections and these are 0.93% in the chemical potential and 
0.72% in the total energy. Due to the cancellation the net corrections are approximately 
of order 1.2% for both chemical potential and total energy over their corresponding TF 
values. Thus the third term which is logarithmic in nature compensates for the correction 
due to the second term in both chemical potential and energy. 

In the next section we present results of our variational calculation along with the 
numbers obtained from these expressions as a further check for our numerical results. 

4 Results and Discussion 

In this section we first report calculations on Rb atoms confined in an isotropic (Ao = 1) 
trap. The choice of isotropic trap is motivated only for the sake of comparison with the 
existing studies []14] . As pointed out earlier, our variational approach is capable of dealing 



with asymmetric traps with cylindrical symmetry or even fully symmetric traps. For the 
symmetric trap we choose trap frequency u /2ir and scaled s wave scattering length a to 
be 77.78Hz and 4.33 x 10~ 3 , respectively. In Table 1 we show chemical potentials, total 
energy and mean radius < r 2 > for particle numbers ranging from 10 3 to 10 9 for both 
GP and MGP cases. For comparison we also show the corresponding MGP results of 



Ref.||14|in parenthesis. It clearly shows that the difference between GP and MGP results 
grow with the increase in particle number. For example, corresponding to N = 10 4 the 
GP numbers are lower than the MGP values by 0.4%. On the other hand for N = 10 8 
the difference is of the order of 1.4%. This is consistent with our estimation of the 
corrections given by Eq.(P2"D and (|2"B|). Moreover, at this point we also note that numbers 
obtained by us for total energy are systematically lower than the corresponding MGP 
results of Ref . ||14|| . The difference between our results and that of Ref.|]T4| increases with 
the number of particles. For example the difference is of the order 2% at N = 10 7 . A 
similar trend is observed for both chemical potential and mean radius. At this point 
it is worth mentioning that our method, being variational, provides an upper bound to 
the ground state energy. It is expected that the numerical value of the ground state 



Table 1: Results for the ground-state of 87 Rb atoms confined in an isotropic trap with 
u/l/27t = 77.78 Hz. Chemical potentials and energies are in units of hu^ and length is 



in units of a>ho- Numbers in the brackets correspond to the results of reference [1 \ 



N 


GP 


MGP 


Mi 


Ei 


^ T ^"max 


Mi 


E x 


*^. ' -^max 


10 3 


3.05 


2.43 


1.65 


3.05 
(3.06) 


2.43 

(2.43) 


1.66 
(1.66) 


10 4 


6.87 


5.05 


2.44 


6.91 
(6.92) 


5.07 
(5.08) 


2.45 
(2.45) 


10 b 


16.90 


12.13 


3.81 


17.02 
( 17.07) 


12.21 
(12.25) 


3.83 

(3.84) 


10 ti 


42.30 


30.24 


6.02 


42.71 
(42.97) 


30.51 
(30.66) 


6.06 

(6.10) 


10 v 


106.18 


75.86 


9.54 


107.53 
(108.75) 


76.75 

(77.48) 


9.61 

(9.74) 


10 8 


266.70 


190.50 


15.12 


270.42 
(275.89) 


193.13 
(196.45) 


15.23 
(15.45) 


10 9 


669.90 


478.50 


23.96 


675.77 


483.81 


24.00 



energy obtained variationally will be typically higher than those obtained by solving the 
MGP equation directly by standard numerical procedures. It is therefore surprising that 
the results obtained by solving the MGP equation numerically in Ref. [14 gives higher 



value of the ground state energy than that obtained by us. On the other hand, numbers 
obtained by the correlated wave function approach in Ref. ||14] are very close to our results 
for iV < 10 6 . For higher values of N, however, the trend is similar to the MGP results. 

Next we compare results obtained by us with those of Ref. [13| . To do this we consider 
BEC of 87 Rb atoms in a trap with angular fequency of 20Hz and a = 7nm. We do 
calculation for iV = 10 9 as for this value of iV quite a drastic energy change is reported in 



Ref.|p.3|. We find that for above parameters the energy per particle corresponding to the 
GP case is lower than the MGP value by 1.38%. This is also consistent with the analytical 

). The corresponding number reported in Ref. |L3| is 33%. We 



estimate given by Eq.(| 
feel that this is a gross overestimate of the energy correction since the first higher-order 
term in the Eq.(|23"D leads to the correction of around 3% over the TF value. But as 
mentioned before the second higher-order term in the Eq.(^3|) being logarthimic, results 
in a negative correction and for above parameters it is 1.6% of TF value. Consequently 
the net correction in energy per particle due to higher-order terms is only 1.4% over the 
TF number. We point out that this number compares well with the result obtained by 
variational calculation. 

In Table 2 we present the numbers for chemical potential and total energies obtained 



Table 2: Analytical estimates for the chemical potenial and the ground-state energy of 
87 Rb atoms confined in an isotropic trap with uj^_/2tt = 77.78 Hz. Chemical potentials 
and energies are in units of Tluj° ± and length is in units of a^. TF represents results 
within Thomas- Fermi approximation (Eqs.(|i"9|) and (|20D) and MTF corresponds to those 
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and (0)) 


N 


TF 


MTF 


Hi 


Ei 


Hi 


Ei 


10 3 


2.66 


1.90 


2.66 


1.90 


10 4 


6.67 


4.76 


6.70 


4.78 


10 5 


16.75 


11.96 


16.87 


12.04 


10 6 


42.07 


30.05 


42.49 


30.33 


10 7 


105.68 


75.49 


107.05 


76.41 


10 8 


265.46 


189.62 


269.28 


192.38 


10 9 


666.81 


476.29 


673.19 


482.42 



via Eq.(p2|) and (p3|), respectively for several values of N. These numbers serve as an 
additional check for correctness of our numerical results especially for large N. This is 
clearly the case as MGP numbers in Table 1 are quite close to the corresponding numbers 
in Table 2 for large N. Furthermore, comparison of Table 1 and II also shows that our 



variational results are closer to the analytical estimates than those reported in Ref. [1TJ] . 

Next to study the density profile we plot the ground-state wave function along x- 
axis in Fig. 1 corresponding to both MGP (solid line) and GP (dotted line) cases for 
three different values of N: N = 10 7 , 10 8 and 10 9 . It clearly shows that for all N the 
difference between the GP and the MGP wave functions arises at the region where it 
reaches maximum, that is at the bottom of the potential well. Moreover, for N = 10 7 
(Fig. la) and 10 8 (Fig. lb) the value of wave function near origin corresponding to 
the GP case is higher than that of MGP case. The difference is of the order 1% only. 
However, the situation becomes just opposite for the case of N = 10 9 (Fig. lc). The 
reason for such change in wave function profile is the contribution of logarithmic term in 
the interatomic energy becoming appreciable at this N. As is already discussed that this 
term introduces negative correction to the interatomic energy and potential. Thus the 
logarithmic term gives rise to an attractive potential for the atoms in BEC at large N. 
Consequently, inclusion of this term in the interatomic energy (potential) brings atoms 
closer to the bottom of potential leading to the fact that more atoms are being found 
near origin than the GP case. Finally it is important to note here that the tail region of 
the wave function is not affected by the higher-order terms in the interatomic potential. 

After having discussed results for isotropic trap we now describe the results for 
anisotropic trap which is more relevant from experimental point of view. For this cal- 
culation we employ the numbers for the asymmetry parameter and the axial frequency 
corresponding to the experiment of Anderson et al. ||TJ. Accordingly, Aq = v8 and 



Table 3: Results for the ground-state of 87 Rb atoms confined in anisotropic trap with 
Ao = Vo and <jj\j2tx = 77.78 Hz. Chemical potentials and energies are in units of Tvjj° ± 
and length is in units of a,ho- 



N 


GP 


MGP 


/ii 


Ex 


fix 


Ex 


10 3 


4.79 


3.85 


4.80 


3.86 


10 4 


10.53 


7.78 


10.58 


7.82 


10 5 


25.66 


18.47 


25.87 


18.60 


10 6 


64.13 


45.86 


64.84 


46.33 


10 7 


160.95 


114.99 


163.15 


116.48 


10 8 


404.24 


288.75 


409.72 


292.66 


10 9 


1015.38 


725.27 


1016.59 


729.65 



uj®/2tt = 220Hz. The value of s-wave scattering length is same as that of isotropic case 
considered earlier. We present these results in Table 3. In this Table we compare the 
MGP results with the corresponding GP numbers. Similar to the isotropic case here 
also the virial relation Eq.flT4|) is satisfied up to 5-th decimal place. The trend in both 
chemical potential and total energy with increase in the number of particles is similar to 
that of isotropic case. The difference between the MGP and GP numbers are of the same 
order as that of isotropic case. Moreover, in the axially symmetric case also we find that 
our variational numbers are quite close to the numbers obtained via Eq.(|2"2"|) and (p3|) 
especially for N > 10 6 . 

Thus we conclude that the variational results obtained for anisotropic trap are also 
quite accurate. Therefore, we demonstrate that by using variational approach described 
in this paper along with a judicious choice of ansatz for the ground-state wave function 
it is possible to obtain reasonably accurate results for the to properties of BEC, even 
beyond mean-field approach, without much of a computational effort. 



5 Conclusion 

In this paper we have studied the properties of Bose gas confined in both isotropic and 
axially symmetric potential going beyond GP or mean-field approximation by taking 
higher-order terms in the interatomic interaction energy. We have used the variational 
approach to solve the MGP equation for wide range of particle numbers. We have verified 
our results using the generalized virial relation as well as by making analytic estimate 
of the corrections introduced by higher-order terms. These higher-order terms lead to 
correction in the total energy, chemical potential and other physical properties of BEC. 
The magnitudes of these corrections are of the order of 1% — 2% even for very large 
N. However, there is qualitative change in the density profile of the condensate due to 
presence of the logarthimic term in the interaction energy. We also critically examine the 
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results in the literature and compare them with our numbers. Here we emphasize that 
the variational method employed by us gives quite accurate results with considerable 
computational ease. 

It is well known that, corrections of this order are difficult to detect in the exper- 
imental H situation. However, these small changes are also reflected in the collective 
excitation frequencies of BEC and these quantities can be measured with greater ac- 
curacy. Motivated by this we are now applying the variational approach and sum rule 
method of response theory of many-body systems |2(J |21| to calculate collective excitation 



frequencies. These results will be presented in our future publication. 
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Figure 1: Plot of the wavefunction along transverse direction for different values of N: 
(a) N = 10 7 , (b) N = 10 8 and (c) N = 10 9 . The solid and dashed curves correspond 
to solutions of GP and MGP equations respectively. The numbers along Y-axis are 
multiplied with a scaling factor of 1000 
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